Baysian Data Analysis: Final Project

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

library(rstan)
library(bayesplot)
library(sp)
data(meuse)
library('plotly')
library(lattice)
source("~/Documents/BU_2019_Fall/bslm.R")

compute_distances <- function(x,y){
  n<- length(x)
  D = matrix(0,nrow=n,ncol=n)
  for(i in 1:n){
    for(j in 1:n){
        p1 = c(x[i],y[i])
        p2 = c(x[j],y[j])
        D[i,j] <- sqrt(sum((p1 - p2)^2))
    }
  }
  return(D)
}

D <- compute_distances(meuse$x,meuse$y)
p <- plot_ly(z=~D)
layout(p,title = "Distances Between Sites", xaxis=list(title = 'Site'),yaxis=list(title = 'Site'))

No trace type specified:
  Based on info supplied, a 'heatmap' trace seems appropriate.
  Read more about this trace type -> https://plot.ly/r/reference/#heatmap
No trace type specified:
  Based on info supplied, a 'heatmap' trace seems appropriate.
  Read more about this trace type -> https://plot.ly/r/reference/#heatmap

p<-plot_ly(x= meuse$x,y=meuse$y,z = ~log(meuse$lead))
layout(p,title = "Log lead concentration as a function of location", scene = list(xaxis=list(title = 'x location'),yaxis=list(title = 'y location'),zaxis=list(title = 'Log lead')))

No trace type specified:
  Based on info supplied, a 'scatter3d' trace seems appropriate.
  Read more about this trace type -> https://plot.ly/r/reference/#scatter3d
No scatter3d mode specifed:
  Setting the mode to markers
  Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
No trace type specified:
  Based on info supplied, a 'scatter3d' trace seems appropriate.
  Read more about this trace type -> https://plot.ly/r/reference/#scatter3d
No scatter3d mode specifed:
  Setting the mode to markers
  Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode

Fix tau2 to be 1 to understand impact of phi.

phi=1
H <- 1*exp(-D/phi)
p <- plot_ly(z=~H)
layout(p,title = "Spatial Correlation for Phi = 1", xaxis=list(title = 'Site'),yaxis=list(title = 'Site'))

No trace type specified:
  Based on info supplied, a 'heatmap' trace seems appropriate.
  Read more about this trace type -> https://plot.ly/r/reference/#heatmap
No trace type specified:
  Based on info supplied, a 'heatmap' trace seems appropriate.
  Read more about this trace type -> https://plot.ly/r/reference/#heatmap

phi=300
H <- 1*exp(D/phi)
p <- plot_ly(z=~H)
layout(p,title = "Spatial Correlation for Phi = 300", xaxis=list(title = 'Site'),yaxis=list(title = 'Site'))

No trace type specified:
  Based on info supplied, a 'heatmap' trace seems appropriate.
  Read more about this trace type -> https://plot.ly/r/reference/#heatmap
No trace type specified:
  Based on info supplied, a 'heatmap' trace seems appropriate.
  Read more about this trace type -> https://plot.ly/r/reference/#heatmap

phi=3000
H <- 1*exp(-D/phi)
p <- plot_ly(z=~H)
layout(p,title = "Spatial Correlation for Phi = 3000", xaxis=list(title = 'Site'),yaxis=list(title = 'Site'))

No trace type specified:
  Based on info supplied, a 'heatmap' trace seems appropriate.
  Read more about this trace type -> https://plot.ly/r/reference/#heatmap
No trace type specified:
  Based on info supplied, a 'heatmap' trace seems appropriate.
  Read more about this trace type -> https://plot.ly/r/reference/#heatmap

xyplot(log(lead) ~ dist.m | soil, data = meuse,
       panel = function(...) {
         panel.xyplot(type = "p", ...)
         panel.lmline(lty = 2, ...)
       })

xyplot(log(lead) ~ dist.m , data = meuse,
       panel = function(...) {
         panel.xyplot(type = "p", ...)
         panel.lmline(lty = 2, ...)
       })

Set up design matrices

y <- log(meuse$lead)
X <- model.matrix(~ -1 + soil + soil:dist.m, data = meuse)

J <- length(levels(meuse$soil))
n <- nrow(X)
p <- ncol(X) / J
site_j <- aggregate(. ~ soil, data = meuse, length)
soil_site <- data.frame(group = site_j$soil,
                             coef = c(rep("_intercept", J), rep("_distance", J)))
X_beta <- model.matrix(~ -1 + coef , data = soil_site)

Metropolis-Gibbs sampling procedure

nchains <- 4
nsims <- 5000
warmup <-1000
gibbs_chain <- matrix(0,ncol = 12,nrow = nsims - warmup)
param_names <-  c("BetaA0","BetaB0", "BetaC0","BetaA1","BetaB1","BetaC1",
                           "Alpha1", "Alpha0", "Kappa2", "Sigma2", "Tau2", "Phi")
colnames(gibbs_chain)<- param_names
nparams <- length(param_names)
sims <- mcmc_array(nsims - warmup, nchains, param_names)

# Initialize starting values for chain using randomly sampled measurements from data and the corresponding posterior estimates for sigma

# Initialize starting values for alpha as MLE estimates for log(lead) ~ f(distance to river), not controlling for soil type.  

n        <- length(y)
beta_fit <- lm(y ~ X - 1)
beta     <- coef(beta_fit)

# Assuming equal weight of deviance for fit of log lead concentration on error in estimation and error from spatial component.
sigma2   <- deviance(beta_fit) / (2*n)
tau2     <- deviance(beta_fit) / (2*n)

alpha_tau <- 10
beta_tau  <- 2

alpha_kappa <- 10
beta_kappa  <- 4

alpha_fit <- lm(beta ~ X_beta - 1)
kappa2    <- deviance(alpha_fit) / nrow(X_beta)
alpha_phi <- 10
phi       <- 1/rgamma(1,alpha_phi,alpha_phi)

# Flat prior on alpha
alpha_prior <- list(mean=rep(0,2),precision=rep(0,2))

H <- function(phi){
  return(exp(-D/phi))
}

phi_jumping <- function(phi_mn){
  phi <- rnorm(1,mean=phi_mn, sd=0.1)
  return(phi)
}

logposterior  <-  function(phi,sigma2,delta,beta,kappa2,tau2,alpha) {
  logpost <-0
  logposterior<- logpost - log(sqrt(sigma2^n)) + (-1/(2*sigma2))*crossprod(y-X %*% beta-delta) - log(sqrt(det(tau2*(H(phi))))) - (1/2) * t(delta) %*% solve(tau2 * H(phi)) %*% delta -
    log(sigma2)  - (alpha_phi+1)*log(phi) - (alpha_phi/phi) - log(sqrt(kappa2^6)) - (1/(2*kappa2)) * crossprod(beta- X_beta %*% alpha) 
   - (alpha_tau+1)*log(tau2) - (beta_tau/tau2) - (alpha_kappa+1)*log(kappa2) - (beta_kappa/kappa2)
       
  return(logposterior)
}

for (chain in 1:nchains) {
  gibbs_chain <- matrix(0,ncol = 12,nrow = nsims - warmup)
  for (it in 1:nsims) {
    # 1. Sample from posterior on alpha
    alpha <- bslm_sample1(beta, X_beta, sqrt(kappa2), alpha_prior)
    
    # 2. Sample from delta
    delta <- rnorm(n,mean=rep(0,n),sd = chol((tau2) * H(phi)) )
    
    # 3. Beta
    beta  <- bslm_sample1(y - delta, X, sqrt(sigma2),
                         list(mean = X_beta %*% alpha, precision = 1 / kappa2))

    # 4. Sigma2
    sigma2  <- (n * (1/n * crossprod(y - X %*% beta - delta))/rchisq(1,n))
    sigma2  <- as.vector(sigma2)
    # 5. Tau2
    nu <- alpha_tau + n/2
    s2 <- (2*beta_tau + t(delta) %*% solve(H(phi)) %*% delta)/(n+2*alpha_tau)
    tau2    <- (nu * (s2/rchisq(1,nu)))
    tau2    <- as.vector(tau2)
  #  tau2 <- 10
    # 6. Kappa2
    nAlpha <-6
    nu <- alpha_kappa + nAlpha/2
    s2 <- (2*beta_kappa+ crossprod(beta - X_beta %*% alpha))/(nAlpha+2*alpha_kappa)
    kappa2 <- (nu * (s2/rchisq(1,nu)))
    kappa2 <- as.vector(kappa2)
    
    # 7. Phi Random walk
    phi_star <- phi_jumping(phi)
    
    logR <- logposterior(phi_star,sigma2,delta,beta,kappa2,tau2,alpha) -          logposterior(phi,sigma2,delta,beta,kappa2,tau2,alpha)
     if (logR >= 0 || log(runif(1)) <= logR) {
          phi <- phi_star
     }
    
    if(it > warmup){
      sims[it - warmup, chain, ] <- c(beta,alpha,kappa2,sigma2,tau2,phi)
    }
    
  }
}

monitor(sims)
mcmc_trace(sims)
mcmc_acf(sims)
mcmc_dens_overlay(sims)
t(apply(sims[,1,,],2, function(x) quantile(x, c(.025,.25,.5,.75,.975))))

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.